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Abstract 

In this article we introduce an asymptotic preserving scheme designed to compute 
the solution of a two dimensional elliptic equation presenting large anisotropies. We 
focus on an anisotropy aligned with one direction, the dominant part of the elliptic 
operator being supplemented with Neumann boundary conditions. A new scheme is 
introduced which allows an accurate resolution of this elliptic equation for an arbitrary 
anisotropy ratio. 



1 Introduction 

The objective of this paper is to introduce an efficient and accurate numerical scheme to 
solve a strongly anisotropic elliptic problem of the form 

r -V • (AV0) = / , inn 

(1) 

[ (/) = on dilo , dz<l> = on dfl^ , 

where C or C M"^ is a domain, with boundary dil ~ dilu U dflz and the diffusion 
matrix A is given by 

■ A± 

The terms A±_ and Az are of the same order of magnitude, whereas the parameter < e < 1 
can be very small, provoking thus the high anisotropy of the problem. In the present paper 
the considered anisotropy direction is fixed and is aligned with the z-axis of a Cartesian 
coordinate system. The method presented here is extended in some forthcoming works to 
more general anisotropies [S]. 

Anisotropic problems are common in mathematical modeling and numerical simulation. 
Indeed they occur in several fields of applications such as flows in porous media [3ul6j , semi- 
conductor modelling [53], quasi- neutral plasma simulations [TT], image processing P^I28| . 
atmospheric or oceanic flows [27], and so on, the list being not exhaustive. More specifically 
high anisotropy aligned with one direction may occur in shell problems or simulation in 
stretched media. The initial motivation for the present work is closely related to magnetized 
plasma simulations such as atmospheric [18l[2T] or inertial fusion plasmas [7l[12] or plasma 
thrusters [I]. In this context, the medium is structured by the magnetic field. Indeed, the 
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motion of charged particles in planes perpendicular to the magnetic field is governed by a 
fast gyration around the magnetic field lines. This explains the large number of collisions 
the particles encounter in the perpendicular plane, whereas the dynamic in the parallel di- 
rection is rather undisturbed. As a consequence the particle mobilities in the perpendicular 
and parallel directions differ by many orders of magnitude. In the context of ionospheric 
plasma modelling [SlITT] , the ratio of the aligned and transverse mobilities (denoted in this 
paper by e~^) can be as huge as ten to the power ten. The relevant boundary conditions in 
many fields of application are periodic (for instance in simulations of tokamak plasmas on a 
torus) or Neumann boundary conditions (see for instance [5] for atmospheric plasmas). The 
system ([T]) is thus a good model to elaborate a robust numerical method. 

The main difficulties with the resolution of problem ([T]) are of numerical nature, as solving 
this singular perturbation problem for small < e ^ 1 is rather delicate. Indeed, replacing 
in the anisotropic elliptic equation e by zero, yields an ill-posed problem, which has an 
infinite number of solutions (namely all functions which are constant in the z-direction) . 
This feature is translated in the discrete case (after the discretization of the problem) into a 
linear system which is very ill-conditioned for e <C 1, due to the different order of magnitudes 
of the various terms. As a consequence standard numerical methods for the resolution of 
linear systems lead to important numerical costs and unacceptable numerical errors. 

More generally, this numerical difficulty arises when the boundary conditions supplied to 
the dominant 0(l/e) operator lead to an ill-posed problem with a multiplicity of solutions. 
This is the case for Neumann boundary conditions, but also of periodic boundary conditions. 
If instead, the boundary conditions are such that the dominant operator gives a well-posed 
problem with a unique solution, this difficulty vanishes as the leading operator alone will 
suffice to completely determine the limit solution. In this case, one can resort to standard 
methods. This is the case of Dirichlet or Robin boundary conditions. In spite of the fact that 
the problem addressed in the present paper arises only with specific boundary conditions, it 
has a considerable impact in many physics problem, such as plasmas, geophysical flows, plate 
and shells, etc. In this paper, we will focus on Neumann boundary conditions because they 
represent a larger range of physical applications, but we could address periodic boundary 
conditions in a similar way. 

Numerical methods for anisotropic elliptic problems have been extensively investigated 
in the literature. Depending on the underlying physics, distinct numerical methods are 
developed. For example domain decomposition (Schur complement) and multigrid tech- 
niques, using multiple coarse grid corrections are adapted to anisotropic equations in [T4j[22] 
and [T31I25] . For anisotropy aligned with one (or two directions) , point (or plane) smoothers 
are shown to be very efficient [23] . A problem very similar to ([T]) is addressed in [15] , treated 
via a parametrisation technique, and seems to give good results for rather large anisotropy 
ratios. However, these techniques are only developed in the context of an elliptic operator 
with a dominant part supplemented with Dirichlet boundary conditions. 
An alternative approach for dealing with highly anisotropic problems is based on a mathe- 
matical reformulation of the continuous problem, in order to obtain a more harmless prob- 
lem, which can be solved numerically in an uncomplicated manner. In this category can 
be situated for example asymptotic models, describing for small values of the asymptotic 
parameter e the evolution of an approximation (f> of the solution of ([1]) [5l|20]. However, 
these asymptotic models are precise only for £ ^ 1, and cannot be used on the whole range 
of values covered by the physical parameter e. Thus model coupling methods have to be 
employed. In sub-domains where the limit model is no longer valid, the original model has 
to be used; which means that a model coupling strategy has to be developed. However the 
coupling strategy requires the existence of an area where both models are valid and still 
demands an accurate numerical method for the resolution of the original model (i.e. the 
anisotropic elliptic problem) with large anisotropics. This can be rather undesirable. 
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In this paper, we present an original numerical algorithm belonging to the second approach. 
A reformulation of the continuous problem ([1]) will permit us to solve this problem in an 
inexpensive way and accurately enough, independently of the parameter e. This scheme 
is related to the Asymptotic Preserving numerical method introduced in |19j . These tech- 
niques are designed to provide computations in various regimes without any restriction on 
the discretization meshes and with the additional property to converge towards the solu- 
tion of the limit problem when the asymptotic parameter goes to zero. The derivation of 
such Asymptotic Preserving methods requires first the identification of the limit model. For 
singular perturbation problems, a reformulation of the problem is required in order to de- 
rive a set of equations containing both the initial and the limit model with a continuous 
transition from one regime to another, according to the values of the parameter e. This 
reformulated system of equations sets the foundation of the AP-scheme. Other singular 
perturbations have already been explored in previous studies, for instance quasi-neutral or 
gyro-fluid limits [101112]. These techniques have been first introduced for non-stationary 
systems of equations, for which the time discretization must be studied with care in order 
to guarantee the asymptotic preserving property. For the anisotropic elliptic equation in- 
vestigated in this article, we only need to precise the reformulated system and provide a 
discretization of this one. 

The outline of this paper is the following. Section [2] of this article presents first the ini- 
tial anisotropic elliptic model. In the remainder of this paper, it will be referred to as 
the Smgiilai- Perturbation model (P-model). The reformulated system (referred to as the 
Asymptotic Preserving formulation or AP-formulation) is then derived. It relates on a de- 
composition of the solution (/)(a;, z) according to its mean part </>(a;) along the z coordinate 
and a fluctuation (f>'{x, z) consisting of a correction to the mean part needed to recover the 
full solution. The mean part (j){x) is solution of an e-independent elliptic problem, and the 
fluctuation (j)'(x,z) = (f>{x,z) — 4'{x) is given by a well-posed e-depcndent elliptic problem. 
The advantage is that the e-dependent problem for the fluctuation is well-posed and solvable 
in an inexpensive way, and this uniformly in e. In the limit e the AP-formulation re- 
duces to the so called Limit model (L- model) , whose solution is an acceptable approximation 
of the P-model solution for £ ^ 1. The present derivation is carried out in the framework 
of an anisotropy aligned along one axis of a Cartesian coordinate system. In the context of 
magnetized plasma simulations, this initial work is extended in a forthcoming work for the 
three dimensional case in curvilinear coordinates, designed to fit a more complex magnetic 
field topology (i.e. anisotropy direction) [6j. The main constraints of this method reside in 
the construction of the mean part which necessitates the integration of the solution along 
the anisotropy direction. This operation is easily carried out in the context of coordinates 
adapted with the anisotropy direction. However, an extension of the techniques presented 
here is currently developed for non-adapted coordinates [9]. 

Section [3] is devoted to the numerical implementation of the AP-formulation. Numerical re- 
sults are then presented for a test case, and the three approaches (AP-formulation, straight 
discretization and resolution of the P-model and L-model) are compared according to the 
precision of the approximation for different values of e. In section |4] we shall rigorously anal- 
yse the convergence of the AP-scheme. Error estimates will be established which underline 
the advantages of the AP-scheme as compared to the initial Singular Perturbation model 
and the Limit model. 

Current research directions arc concerned with the adaptation of the present technique 
to the case of arbitrary spatially varying anisotropics, without adaptation of the coordinate 
system to the direction of the anisotropy. These developments will allow the treatment on 
nonlinear problems, when the diffusion tensor (and its principal directions) depend on the 
solution itself. This treatment will involve iterative methods which, at each iterate, will 
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reduce the problem to the sohition of a Unear anisotropic diffusion problem. 



2 The asymptotic preserving formulation 

For simplicity we shall consider in this paper the two-dimensional problem, posed on a 
rectangular domain O = fi^; x fi^, where il^; C M and fi^ C M are intervals. The ideas 
exposed here can be extended without any problems to the more physical three-dimensional 
domain, with two transverse directions (x, y) and an anisotropy direction aligned with the 
z-direction. In this section we introduce the Singular Perturbation Model, the Limit Model 
and the Asymptotic Preserving formulation. 



2.1 The Singular Perturbation Model (P-model) 

The main concern of this paper is the numerical resolution of the following anisotropic, 
elliptic problem, called in the sequel Singular Perturbation Model 

r -V-(AV0) = /, in n, 

— — = on ilx y. dilz , = on d^lx x fi^ . 
I. dz 

The anisotropy of the media is modeled via the definition of the diffusion matrix A 

"o ,1 ) . 

where A±(x, z) and Az{x, z) are given functions with comparable order of magnitudes. The 
source term /(cc, z) is given and the parameter e is small compared to both A± as well as 
Az- The medium becomes more anisotropic as the value of e goes to zero. 



2.2 The limit regime (L-model) 

In this section we establish that in the limit e the solution of the perturbation model 
converges towards (/), solution of the L-modcl defined by 




f{x) , in n 



X ) 



where overlined quantities designate averages over the z-coordinate : 

fix) = ^J fix,z)dz. 



First we can rewrite the P-model as 



dz 

and integrating along the z-coordinate gives 



(4) 



on fix X dflz 1 = on dflx x ft^ , 
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This equation holds for any e > 0. Now, letting formally e tend to zero in ^ yields the 
reduced model (R-model) 

-1(a.'4)-o, in n, 



(i?) <; ^ (7) 

— = on fix 'X dilz J = on dfl^ x fl^ . 
oz 

The functions verifying this ill-posed R-model are constant along the z-coordinate. Thus 
including this asymptotic limit property into equation ^ gives rise to the L-model (0]), 
verified by the solution of the Singular Perturbation model in the limit e 0. 

Remark 2.1 The L-model is the singular limit of the original P-model ([2]). It provides an 
accurate approximation of the P-solution only for small values of e. The P-model is valid 
for all < e < I, but numerically impracticable for e <^ 1. Indeed working with a finite 
precision, the asymptotic model degenerates into the R-model defined by ^ as e vanishes. 
This R-model is ill-posed since it exhibits an infinite amount of solutions (j) = 4>{x), depending 
only on the variable x. This implies that the discretization matrix derived from the P-model 
is very ill-conditioned for small < e ^ 1. This point is addressed by the numerical 
experiments of section \3.S\ Consequently, in a domain where e varies significantly, a model 
coupling method has to be developed in order to exploit the validity of each model, the P- and 
L-model. This can be rather undesirable. In the next section we shall present an alternative 
approach, which is based on a reformulation of the Singular- Perturbation model providing 
a means of computing an accurate numerical approximation of the solution for all values 
< e < 1. 

Remark 2.2 The asymptotics is totally different in the case of Dirichlet boundary con- 
ditions. In this case, the R-model is well posed, with a unique solution, and there is no 
difficulty anymore. Any standard numerical solution of the P-model will converge to that of 
the R-model. In other words, with Dirichlet boundary conditions, the perturbation becomes 
regular and the limit solution is fully determined by the formal limit system. The situation 
and the difficulty addressed in the present paper require that the R-model be ill-posed. This 
is the case with Neumann boundary conditions (which is the framework chosen here) but 
also with periodic boundary conditions, or any other boundary condition which would result 
in an ill-posed R-model. 



2.3 The Asymptotic Preserving reformulation (AP-formulation) 

In order to circumvent the just described numerical difficulties in handling the Singular 
Perturbation model, we introduce a reformulation, which permits a transition from the 
initial P-model to its singular limit (L-model), as e ^ 0. 

For this, we shall decompose each quantity f{x, z) into its mean value f{x) along the z 
coordinate and a fluctuation part f'{x, z). For simplicity reasons let in the following Qx ■= 
{0,Lx) and fl^ := (0,L^). Then 

fix,z) = f{x)+f\x,z), (8) 

with 

f{x) 7- / f{x,z)dz, 
Note that we have the following properties 

/' = o, Wm^dj/sx, Jg^fg + W, (10) 

af/a. = of'/9,, {df/9x)' = df'/a., {fgY^f'g'-W + .fg' + f'9- (H) 



f'{x,z):^f{x,z)~f{x). (9) 
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Taking now the mean of the eUiptic equation ^ along the z-coordinate, we get thanks 
to (fTO]) and pT|) . an equation for the evolution of the mean part (j){x) 



d 



[API) 



d4> 



dx 



d 



dx 



dx 



in , 



(12) 



on dflr 



Substracting from ^ this mean equation gives rise to the evolution equation for the 
fluctuation part (l)'{x, z) 



oz 



dx 



dx 



{AP2) < 



d_ 

dx 



, di_ 

dx 
A\ 



d(j)' 



= on fir X dfl. 



, m fix 



^'dx^ 

y = on dVlx X , 



(13) 



Thus we have replaced the resolution of the initial Singular Perturbation model ([5|) by the 
resolution of the system (fT2|l -()13 p . which will be done iteratively. Starting from a guess 
function cf)' , equation (|12p gives the mean value 4>{x)^ which inserted in (|13p shall give the 
fluctuation part 4>'{x, z) and so on. 

The constraint (/)' = in ([T3|) (which is automatic for e > 0, as explained in Remark 
12. 3p has the essential consequence that the conditioning of the discretizcd system becomes 
£-independent, because the problem ()13p reduces in the limit e ^ to the system 



dz 



dz 

on ^x 'X dQ 



, in n, 



dz 

0' = in Qx , 



on dQx X flz 



(14) 



which is uniquely solvable, with the solution cp' = 0. Inserting this solution in (|12p . we 
conclude that the solution of the AP formulation converges for e — > towards the mean 
value part (t>{x), computed thanks to the Limit model 



d 



dx 



dcj) 



dx 



in fix ■ 



(15) 



= on dflx . 



The AP reformulation (|T2|) -(|T3ll is equivalent to the Singular Perturbation problem (jSj) 
and is therefore valid for all < e < 1. This new formulation guarantees that, working with 
a finite precision arithmetic, the computed solution converges in the limit e ^ towards 
the solution of the limit model ([4]). This is a huge difference with the original Singular 
Perturbation model which degenerates into an ill-posed problem. Thus, by using the AP- 
formulation, we expect the computation of the numerical solution to be accurate, uniformly 
in s. 

For the detailed mathematical proofs, we refer to the next section. 
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Remark 2.3 The condition 0' = m US]) holds automatically for e > 0, since the right- 
hand side has zero average along the z-coordinate. Indeed, let tp be the solution of 



d 



d 



^{A._^]-e^lA,^]+e^lA\^ 



dz J 



dx 



di' 



dx 



d 



dx 



dx 



£5 



in . 



dip 

— — = on r^a: X dflz , 
oz 



(16) 



'0 = on d^lx X 



with g' = 0. Taking the average along z, we get 

dip' 



d_ 

dx 



A, 



dx 



0, 



in 51-, 



■0 = on dQx , 



and thus ip = 0, which is nothing but the constraint added in I113\) . 

The computations of the fluctuating part (p' via the equation (|13p requires the discretization 
of an integro- differential operator. This means that the discretization matrix will contain 
dense blocks. However, using (|12p the system (AP2) can he rewritten as 



{AP2') 



d_ 

dz 



dcl/_ _ 
dz 



A._ 



d_ 

dx 



A, 



dx 



dx 



A 



d(p_ 
dx 



in 51 . 



(17) 



on ilx dflz 
in ilx . 



on dflx X ^z 



In this expression the right-hand side has no longer zero mean value along the z-coordinate, 
but the integro- differential operator has disappeared. The associated discretization matrix is 
thus sparser than that obtained from the system (fT2|) . Systems il ^ )- (T^) and U^ )-[Tl\l are 
equivalent. 



2.4 Mathematical study of the AP-formulation 

We establish in this section the mathematical framework of the AP-formulation p^ - p^ 
and study its mathematical properties. Let us thus introduce the two Hilbert-spaces 

V := {tP{-, •) G H\n) / V = on 911^ x Q^} , W := {V'(-) e H\nx) / ip = on 911^} , 

with the corresponding scalar-products 

{(p, ip)v ■■= e{dx(p, dxip)L^ + {dz(p, dzip)L'2 , (</>, ip)w ■= {dxCp, dxtp)L'^ , (18) 

and the induced norms || • ||v, respectively || • ||vv- For simplicity reasons, we denote in 
the sequel the scalar-product simply by the bracket (•,•). Defining the following bilinear 
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forms 



"'0 



ai ((/)', ■0') ■—I / A^{x^ z)——{x, z)— — {x,z)dxdz 



"'0 



9a; 9a; 



02(0,'/') := / Ax(a:)^(a;)^(a:)dx, 







9x dx 

c(0',V;) := /"^^ f"^^ A'Ax,z)^ix,z)^{x)dxdz, ^^^^ 



^0 



b{P,ij') ~ / P{x) / il;'{x,z)dzdx, 
Jo Jo 

a{cp', := ao (0', V') + eai (</>', - ed{^' , , 

permits to rewrite the AP system (IT2t - p^ under the weak form 

a2(0>) ;^c(0',Vi) , Vt^gW, 

6(Q, 0') = , VQ e W , 

where (/)'(a:;, z) £ V, (t){x) e W as well as P{x) G W arc the unknowns and V'' G "1^- "0 G W 
and Q G W the test functions. It can be observed that the constraint 0' = was introduced 
via the Lagrange multiplier P. We will see in the next theorem that the weak formulation 
(|20|) is equivalent for e > to the system 

a2(0>) = (/"Vi)-^c(0'>) , V^'GW, (21) 

a (0', yj') = £(/', 7^') - ec {yj',^ , V^' e V , (22) 
where the explicit constraint ((>' = does not appear. Let us assume in the sequel 

Hypothesis A Let the diffusion functions A± E i°°(il) and A^ G L°°(il) satisfy 
< c_L < A±{x, z) < M± , < < A^{x, z) < , fa.a. (x, 2) G , 
wiit/i some positive constants c±, Cz, M±, M^- Let moreover f G L'^{fl). 

The next theorem analyzes the well-posedness of the AP-formulation. 



Theorem 2.4 For every e > the problem Ii21 \) - l2^) admits under Hypothesis A a unique 
solution [4>'^,4>^) G V X W, where (f>e ■= 4''e + 0e ^■s the unique solution of the Singular 
Perturbation model {^i- The function 0^ has zero mean value along the z-coordinate, i.e. 
(j)'^ = for every e > 0. 

Consequently, (4>'^, 4>^) G VxW is the unique solution of i21\) - i2S\) if and only if {<j)'^,(j)^, P^) G 
V X W X W is a solution of the AP-formulation I120\}. In this last case, we have Pg — 0. 
Finally, these solutions satisfy the bounds 



An asymptotic preserving scheme for strongly anisotropic elliptic problems 



9 



with an e-independent constant C > 0. In the limit e — there exist some functions 
(i^q,(/)q) e V X yV, such that we have the following weak convergences in 

and the strong convergences 

<?!>£ ^£-+0 00 L'^i^), 9^01 ->£->o <92 0o in L^iVl), 0^-^.^0 00 *^ L^{9.x), 

where 0q = and 0q is the unique solution of the Limit model 

Proof: The Singular Perturbation model ([5]) and the Limit model ^ are standard elliptic 
problems and posses under Hypothesis A (and for every e > 0) unique solutions 0^ € V, 
respectively € W. It is then a simple consequence of the decomposition ([5]), that the prob- 
lem ([2T|) -(|22 |) admits a unique solution (0^,0^) S V x W, where 0£(a;) '■= Jq"" (f^eix, z)dz 
is the mean and 0^ := 0£ — 0^ the fluctuation part. Thus we have also 0^ = 0. This property 
can also be understood from the fact that the right-hand side of (fT^ , denoted in the sequel 
by 3 

g{x, z) := f'{x, z) + ^)|~(^)^ ' 

has zero mean value along the z-coordinate. Indeed, taking in test functions ip'{x) £ V 
depending only on x, yields immediately that 0^ = for all e > 0. 

Standard stability results for elliptic problems yield now the e-indepcndent estimate for the 
solution of the Singular Perturbation model ^ 

< \\d.A\\l2^n) + -||5.0£lli2(n) < ^^ll/llL(n) , 

implying that 110^11^1(0^ < C'll/lli2(n) and ||0il|^i(jj) < C||/| , with a constant C> 
independent of e > 0. Thus there exist some functions (0o, 0o) G V x W, such that, up to a 
subsequence 01. ^£^o 0o in H^{fl) and 0^ ^£^o 0o in H^(flx). Hence we have 

I I (j)'^{x, z)'tp{x, z)dx dz ^e-*o / / (I)'q{x, z)'tp{x, z)dx dz , V?/' G V . 
Jo Jo Jo Jo 

Taking here "0(2;) G V depending only on the x-coordinate, we observe that the feature 
0'^ = yields the crucial property of the limit solution 0'q = 0. Passing now to the limit 
£ ^ in ((22l) . we get that 0q is solution of 

ao(0[,, V') = 0, VV^'gV, with 0^=0 in il^ , 

which is the weak form of and implies 0q = 0. Finally, passing to the limit in (PT|) . 
yields that 0q is the unique solution of the Limit model (jlj. Because of the uniqueness of 
the limit (0g, 0g), we deduce that the whole sequence (0^, 0^) converges weakly towards this 
limit. To conclude the first part of the proof, we shall show the strong convergences. For 
this, taking in ([22|) 0^ as test function and passing to the limit e ^ 0, yields (9z0g ^ in 
L^(ri). As 0g G V and 0^ = 0, the Poincare inequality 

II0;||l^ <q|5.0;iU2, 

is valid and implies that 0^ ^ in L^{fl). The convergence 0^ 0o in i^(rix) is immediate 
by compacity. It remains finally to prove the equivalence between (|20p and (|2ip - (P^ . This 
is immediate. Indeed, if (0^, 0^) G V x W is solution of pT|) - ((22)) . then (0^, 0^, 0) is solution 
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of (Uni)- And if (0^, 0^, Pe) e V X W X W satisfies then = (obvious by taking as 
test function in ((20|) ip'{x) £ V depending only on x) and solves hence ((2T|) - ((22|) . ■ 

The subject of the next section will be the numerical resolution of the AP-formulation 
(fT2 |) -(fT3 l) (or ((20|) ) and this shall be done iteratively via a fixed-point apphcation. Let us 
thus introduce here the fixed-point map, construct an iterative sequence and analyze its 
convergence. In the rest of this section, the parameter e > shall be considered as fixed. 
Due to the fact that the two systems (|20|) and ([2T|) -([22 |l are equivalent, we shall concentrate 
on the simpler one, i.e. ([2T |) -([22 ll . Let us define the Hilbert space 

U ■.^{i;{;-)eV /4' = 0}, 

associated with the scalar product 

(0,-0)*:= / / A^dz(t)dzijdzdx + e / A^d.j,(j)dx'>pdzdx , 
Jo Jo Jo Jo 

which is equivalent to the scalar product (•, •)v on V, defined by p8|) . 

The fixed-point map T -.U ^ U is defined as follows: With cj)' eU wc associate (p G W, so- 
lution of (pi]) . Then constructing the right-hand side of via this (f> £ W, we define T {(/)') 
as the corresponding solution of (f^^ . Denoting by (0'*, 0*) G V x W the unique solution of 
(|2T|) - (|22)) . we remark by Theorem 12.41 that 0'^ € U and that it is the unique fixed-point of 
the map T. 



Theorem 2.5 Let s > be fixed and let (j)'^ G U be the unique fixed-point of the application 
T : U U constructed as follows 



Then for every starting point 0q G U, the sequence (f)'^. := T{(f>'i^_^) ~ T'^(0g) converges in 
{U, II • II*), and consequently also in {hi, \\ ■ \\v), towards the fixed-point 4>'^ Cz U of T . 



The proof of this theorem is based on the following 



Lemma 2.6 \SI Let {lA,\\-\\<,) be a normed space and T :U ^ U a contractive application, 
i. e. 

\\T{(b)~T{i')\u<U~i,\U, with cjj^i;. 

Then the set of fixed-points of T, denoted by FP{T), is identical with the set of accumulation 
points of the sequences {T'^(0)}fegN, with cj) gU, set which is denoted by AP{T). Moreover, 
these two spaces contain at most one element. 

Proof of theorem 12.51 : 

The linear application T is well-defined. The first step cj)' Gl4 > G W is immediate 

by the Lax-Milgram theorem. For the second step, we remark that for given G W the 
equation 

a(0,^') = £(/',V'')-£c(0',0), V^'GV, (23) 

has a unique solution 6 Cz U. Indeed, we notice first (by taking test functions only depending 
on the x-coordinatc) that 9 — 0. This enables us to consider instead of the variational 
formulation 

mi9,i;') = e{r,i'')^ec{^',^), W-' € V , (24) 
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where the bihnear form a(-, •), which is not coercive, was replaced by the coercive bihnear 
form m(-, •), given by 



eMi 



dxO{x, z)dz 



dxip'ix, z)dz 



dx . 



(25) 



Indeed, due to the property = 0, the two equations (pS)) and are equivalent and this 
time m(-, •) is a continuous, coercive bilinear form, as for all 0' G V we have 



"'0 



A^ld^tp'l dzdx + e 



A^\dx^'\'dzdx > C||V''||v 



Thus the Lax-Milgram theorem implies the existence and uniqueness of a solution 9 Gt( oi 
the continuous problem (|24[) and hence also of problem We have shown by this that 

T is a well-defined mapping. 

Furthermore we know that T admits, for fixed e > 0, a unique fixed-point, denoted by 
(f)'^ € hi. Let us now suppose that we have shown that T is contractive. Then lemma 12.61 
implies that FP{T) = AP{T) = {4>'^}. Thus choosing an arbitrary starting point cf)Q £ 14, 
and constructing the sequence (p'l, := r((/>'j,_x) = T''{(I)q), we deduce that this sequence has 
a unique accumulation point (f>'^ in U. This means that the sequence {(j}'f.}ki£N converges in 
{U, II • ||,) towards 0'^. Due to the fact that || • ||* and || • ||v are equivalent norms, we have 
also the convergence in {U,\ \ ■ \ 

It remains to show that T is contractive. For this let (/)[ ,<j>2 E U he two given, distinct 
functions. Denoting by 0' := (p'l — (p^, '■— 4>i ~ 4>2 (where (pi eW are the corresponding 
solutions of (HIl)) and 9' := T(0i) - T(0y, we have to show that ||6''||* < ||(^'||*. First we 
observe that <f) solves 

a2(^>) = -^c(^',V5), V^eW, 



and 9' is solution of 

a{9',tp') = -ec{ip','4>) , VV-' e V. 
Taking in as test function, gives rise to 



(26) 
(27) 



A^\dx<p{x)\'' dx = - 



1 



Jo 



— / A'^dxcl)'{x,z)dz 



— / A^dx(j)'{x,z)dz 



dx4>{x) dx 
dx4>(x) dx 



< 



Thus 



A±\dxcj){x)\'dx 











' r^A 




z 


Jo Jo 


















Jo Jo 



- 1/2 






1/2 






\dxp\'dx 






Jo 










1/2 


L\dx4 


'\^dzdx 





Equally, taking in ((27)) 9' as test function gives rise to 
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A,\d,e'\^dzdx + e / A±\d^9'\^dz < ~e / A±d^'^ d^O' dzdx 
Jo Jo Jo Jo Jo 



< e 



Jo 





1/2 




1/2 


A±^\dx4>'\^dzdx 












Jo Jo 





(28) 



This last inequality yields 
A^\d^e'\'^dzdx 



A^\d.,<j>\'dx 



- 1/2 




1/2 




/ / A±\d.J'\'^dzdx 






Jo Jo 





el / A^\dxd'\^dzdx < eL^ A±\dx'^\^dx 
'o Jo Jo 

< e A^\dx(l)'\'^dzdx 

Jo Jo 

< /" Y 'Az\d,(j)'\^dzdx + e [ Y 'A±\dx(P'fdzdx. 
Jo Jo Jo Jo 



In this last step we would have the "equality" if and only if J^" J^' Az\dz(j)'\'^dzdx = 0. 
This is however only possible for functions depending exclusively on the x-coordinate, 4>' {x), 
which is in contradiction with the fact that = and (/)' 7^ 0. Thus we have shown that 
||r(0')||* < ||0'||* for 0' 7^ 0, 0' G which means that T is a contractive application on 



3 Numerical discretization and simulation results 

This part of the paper is concerned with the numerical discretization of the AP-schcme 
(|12p - p3p and the comparison of the simulation results with those obtained via the Singular 
Perturbation model ^ and the Limit model (|4|). 

3.1 Discretization 

The numerical resolution of the Asymptotic Preserving system (jl2p -(|13 p is done by means 
of the standard finite clement method. 

Let us recall the variational formulation of the AP-formulation 

f a2 (0, = (/, ^) - j-c (0', Vi) , Vt^ e W , 

I a (0', V') + h{P, Vy ) = V^) - ec ,'4>) , VV-' £ V , (29) 

I &(g,0') = o, vgGW, 

with the notation of section [51 Here (l)'{x, z) £ V, <f>{x) e W as well as P{x) G W are the 
unknowns and ip' G V, ^ G W and Q G W the test functions. 

The introduction of the Lagrange multiplier P{x) was explained in a simplistic manner in 
the preceding sections and will be analyzed in more details in section |31 Due to the equiva- 
lence of ((29)) and (|2T|) - ((22)) . one can comment that the introduction of P{x) is superfluous, 
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but this is not the case for the discretized equations. The property 0' = is indeed au- 
tomatically fulfilled since the right-hand side of equation ((22|) has a zero mean value along 
the z-coordinate. However the discrete implementation of this quantity introduces round-off 
errors which probably will destroy the zero mean value property and justify the introduction 
of the Lagrange multiplier. 

For simplicity reasons wc omitted here the e-indcx of the solution (cj)'^, (p^), the parameter 
e > being considered as fixed. 

To discretize now the system ()29|) we introduce the grid 

= Xo < • • • < < • • • < XN^ + i ^ Lx , = Zi < • • • < Zfc < • • • < Zjv, = 

and denote the cells by /„ := [a;„,x„+i] and Jk ■= [zk,Zk+i]- The finite dimensional spaces 
Vft, C V and Wh C W are constructed as usual, by means of the hat functions (Qi finite 
elements) 



J X ^ Ifi—l ^ 



Xji Xji—i 

Xnix):^{ Xn+i-x ^ ^^^^^ ^ Kk{x):={ 

Xn+l Xji 

, else 



Zk-l ^ J 

z G Jfc_i 



Zk - Zk-l 
Zk+1 - z 



, z € Jk, 



Zk+i — Zk 
, else 



Thus we are searching for approximations cp'/^ E Vh, 4>h G VV/i and Pk € Wh, which can be 
written under the form 

'i>'h{x,z) =^^ankX7i{x)K.k{z) , (j)h{x) ^^PnXnix) , Ph{x) =^-1nXn{x) . 
n— 1 k—1 n—1 n—1 

Inserting these decompositions in the variational formulation (|29p and taking as test func- 
tions the hat-functions Xn and Kk gives rise to the following linear system to be solved in 
order to get the unknown coefficients a„fc, /3„ and 7„ 

(30) 

AQ + e{Ai-D) B \ f a \ _ f v 



where the matrices A2 € R^-x^-, Ao,Ai,D G R^=^^-x^-^- and B G R^-^-x^- correspond 
to the bilinear forms ([TO]) and the right-hand sides are defined by 

W„ := {f,Xn) - 7-0(0^, XrO , V„fc := {f',XnKk) " c(x„Kfc,(^/i) = {9,XnKk) , 

for all n = 1 , • • • , Nx ; = 1 , ■ ' ■ -^z and 

5(x,z):=/'(a;,z) + ^(^Al(x,z)^(.T)^ . (32) 

Solving iteratively the linear systems (j30p -pi p permits finally to get the unknown function 
<ph(x,z) — (ph(x)-\-(p'yJyX,z'). The convergence of the iterations was proved for the continuous 
case in theorem 12.51 and can be identically adapted for the discrete case. 



3.2 Numerical results 



In this section we shall compare the numerical results obtained by the discretization of the 
Singular Perturbation model, the Limit model and the just presented Asymptotic Preserving 



An asymptotic preserving scheme for strongly anisotropic elliptic problems 



14 



reformulation. With this aim, we consider a test case where the exact solution is known. 
Let thus 

(Pf.[x, z) := sm \^—x j + e cos [—z j sm [^—x j , (33) 

be the exact solution of problem where we choose A±{x, y) = ci + xz^ and Az{x, z) — 
Ci + xz, with two constants c\ > 0, C2 > 0. The numerical experiments are performed with 
Lx = Lz = 10 and ci = C2 = Lz- The exact right-hand side / is computed by inserting ([55]) 
in ([5]). We denote by (/)p, c/)]^ and (j>A, respectively, the numerical solutions of the Singular 
Perturbation model ([5]), the Limit model ([U and the Asymptotic Preserving formulation 
p^ - (fT5)) . The comparison will be done in the ^^-norm, that means 

\ 1/2 

11413 - (l)nufn\\2 ^ —7= [ y '\4'e{Xi) - 4>nu7nA'^ ] , (34) 



where 4>num stands for one of the numerical solutions and (j)e{Xi) is the exact solution 
evaluated in the grid point Xi. The index i covers all possible grid indices, reassembled in 
the set G, and is the total number of grid points. The linear systems obtained after the 
discretization of either the P-model, the L-model or the AP-formulation are solved thanks 
to the same numerical algorithm (MUMPS [2]). The purpose here is not to design a specific 
prcconditioner for the resolution of these linear systems, but to point out the efhcicncy of 
the presently introduced AP-mcthod to deal with a large range of anisotropy ratios. 
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(a) Grid with 50 X 50 points. (b) Grid witli 500 X 500 points. 

Figure 1: Absolute error in the Z^-norm between the computed solutions 4>p, ^l^'I'a E^nd 
the exact solution 0e, as a function of e and on different grids. Dashed lines : (S) 
Standard scheme : discretization of the P-model; Stars : (AP) AP-scheme; 
Circles : (L) discretization of the L-modcl. 



As can be seen from Table [T] and Figure [TJ the finite element resolution of the Singular 
Perturbation model is precise only for large < e < 1, whereas the Limit model is accurate 
for small e ^ 1. The range of £- values in which both the Singular Perturbation and the 
Limit models provide an accurate approximation of the solution shrinks as the mesh size is 



refined. For a coarse grid (with 50 x 50 points see figure 1(a) ) this domain ranges from 10^^^ 
to 10~^ while it is reduced to 10~^ — 10~^ for the refined 500 x 500 grid (figure 1(b) ). This 
question is determinant for the development of a model coupling strategy. Indeed it requires 
an intermediate area where both discretizcd models furnish an accurate approximation and 
we observe that for refined meshes this area may not exist. This reduction of the validity 
domain can be explained for both the L-model and P-modcl but for quite different reasons. 
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£ 


10 


1 


10-^ 


10-i 


10-^4 


io-i« 


AP-schemc 


3.4-10"^ 


7.8 • 10"^ 


3.8 • 10-^ 


2.7-10"^ 


2.7-10"^ 


2.7-10-^ 


S-schcmc 


2.8 • 10-^ 


4.5 • 10-^ 


2.8 • 10-^ 


2.7-10-^ 


6.6 • 10-^ 


1.2 


L-modcl 


9.9 


1.0-10^ 


1.0-10-^ 


2.8 • 10-^ 


2.7-10-^ 


2.7 • IQ--^ 



Table 1: Absolute error in the Z°°-norm for the approximation computed thanks to the 
AP-scheme, discretized Singular Perturbation and Limit models (S-scheme and 
L-model) as compared to the exact solution. 



The numerical approximation computed via the Limit model is altered by both the 
discretization error of the numerical scheme and the approximation error introduced by the 
reduction of the initial Singular Perturbation problem to the Limit problem. For coarse 
grids, the global error is rapidly dominated by the scheme discretization error, but as the 
mesh is refined, the approximation error becomes preponderant, as the Limit model is precise 
only for small e-values. The schemes implemented here are of second order, thus when the 
mesh size is divided by ten, the discretization error is reduced by one hundred. The global 
error for the L-model displayed in figure 1(a) does not depend on e as soon as e < 10"'^. 
Below this limit the L-model is able to furnish a better approximation of the solution with 
vanishing e, however the numerical scheme is not precise enough and consequently the global 
error does not decrease. For the refined mesh, this discretization error is lowered by two 
order of magnitudes and the global error is a function of £ as long as its value is greater 
than lO-'^ (Fig. [T(b)| ). 

The analysis for the Singular Perturbation model is quite complementary. The accuracy 
of the approximation provided by the P-model is good for large e-values and deteriorates 
rapidly for small ones. This can be explained by the conditioning of the linear system 
obtained by the P-model discretization. An estimate of the condition number for the matrix 
is displayed in figure [2] for two different grid sizes. This conditioning deteriorates with 
vanishing e-parameter, which is coherent with the fact that, working with a finite-precision 
arithmetic, the Singular Perturbation model degenerates into an ill-posed problem. This 
also explains the blow up of the error displayed in figure [1] as soon as the conditioning of 
the matrix approaches the critical value of the double precision (materialized by the level 
10^^ in Fig. [2]). This limit is reached on more refined meshes for larger e-values (e 10^^^ 
on a 50 X 50 grid and e w 10~^° on a 200 x 200 grid). As expected, the P-model, though 
valid for all e-values, cannot be exploited numerically for small e. The e-region where both 
the P-model and the L-model are accurate all-together, shrinks dramatically with the size 
of the mesh, fact which motivates the development of the AP-method. 
The condition number estimate of the linear system providing the approximation of the 
solution for the AP-scheme is also plotted in Figure [2] The conditioning of the system is 
rather e independent and this is due to the introduction of the Lagrange multiplier, which 
forces the system in the limit to remain well-posed. The accuracy of the AP-schemc is totally 
comparable to the P-model for the large values of e and to the L-model for the smallest ones. 
The AP-formulation is a good tool for computing an approximation for the solution which 
is accurate uniformly in < e < 1 and is therefore of great practical interest. Note that this 
approximation is obtained thanks to an iterative sequence {05;}feeNi constructed with the 
fixed-point mapping T defined in theorem 12.51 The convergence of this iterative process is 
analysed in figure [3] on a 200 x 200 grid for a large value of e. The /^-absolute error between 
the mean respectively the fluctuating parts of the exact solution and the approximation 
provided by the AP-scheme are plotted as a function of the iteration number. The sequence 
is initiated with the zero function. With the iterative process, both components converge 
towards the solution until the precision of the schemes is reached. At this point, after 
roughly 27 iterations, the approximation can not be improved and a plateau is observed. 
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S-scheme (50x50) 
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Figure 2: Condition number estimate for 
the discretization matrices of 
the Standard (S) and AP 
schemes (computed by 
LAPACK [4]) as a function of 
e. Different grids of 50x50 

and 200x200 points and 
different e-values are used. 
Dashed/Plain hues : 200 x 200 
/ 50 X 50 grid ; Stars : 
AP-scheme. 
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Figure 3: The P absolute error between 
the exact solution and the 
numerical approximation 
computed with the 
AP-scheme, as a function of 
the iteration number, with 
£ = 10 and a 200 X 200-mesh. 
Dashed line : mean part of the 
solution; Plain line : 
fluctuating part. 



The convergence of this sequence may be improved thanks to classical relaxation techniques. 
Finally we investigate the positivity of the AP-scheme. With this aim the anisotropic elliptic 
problem is solved with a positive source term, in this case an approximation of the Dirac 
(5-function. This function denoted (5j has a support included in a subset ([—a, a] x [—a, a], 
with < a < 1) of the simulation domainr [—1,1] x [—1, 1]. Two different parameters a are 
chosen, a = 10~^ and a = 10~^. 

The simulation domain is discretized by a 500 x 500 mesh. For the smallest value of 
a the support of the function is reduced to 5 cells in each direction. The source term (5^ 
is normalized, such that the maximal value of 5^ grows with vanishing a-parametcr. In 
table [2] the maxima and minima of the numerical approximations computed by the AP- 
schemc {(f>A) and the discretized Singular Perturbation model {(j)p) are gathered for the 
two source functions (5^. Only large e-values are considered to verify the positivity of the 
numerical approximations. Indeed for very small e the solution is reduced to its mean part 
which is the solution of a classical elliptic problem preserving the maximum principle. This 
means that the relevant question is related to configurations where the fluctuating part 0' 
has a significant contribution to the elliptic problem solution. In this range of large and 
intermediate e values, both approximations are comparable. Only slight differences can be 
observed on the maxima for the smallest e-parameters. The results of table [2] demonstrate 
the positivity of the approximations computed by either the AP-scheme or the Singular 
Perturbation model. 
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e 


102 


10 


1 


10-1 


10-2 


10-3 



7 

O 
1 — 1 


niax(0p) 


77.58 


3.82 


1.63 


8.93 


7.22 


6.93 


niax(0^) 


77.58 


3.82 


1.63 


8.93 


6.89 


6.89 


min(0p) 


1.9 lO-'^ 


2.510-"^ 


2.410-2 


2.410-2 


2.810-2 


2.810-2 


min(0A) 


1.910-^ 


2.510-^ 


2.410-2 


2.410-2 


2.810-2 


2.810-2 




O 
1 — 1 

It 

e 


niax(0p) 


1.8102 


7.1101 


2.6101 


1.2101 


8.29 


7.34 


niax(0^) 


1.8102 


7.1101 


2.6101 


1.2101 


7.14 


7.11 


min(0p) 


1.610-^ 


2.510-3 


2.410-2 


2.810-2 


2.810-2 


2.810-2 


min((?!)^) 


1.6 lO-'^ 


2.510-3 


2.410-2 


2.810-2 


2.810-2 


2.810-2 



Table 2: Maxima and minima of the numerical solutions computed thanks to the 

AP-scheme {(Pa) and the Singular Perturbation model (0p). The elliptic problem 
is solved with the Dirac function as a source term on a 500 x 500 mesh. 



4 Numerical analysis of the AP-scheme 

In this last part of the paper we shall concentrate on the numerical analysis of the Qi finite 
element scheme introduced in section [3. II for solving 



^ ( A ^] - e— ( A ^] + e— [ A' ^ \ - e in O 

dz \ ^ dz ) ^ dx \ ^ dx ) ^ dx \ dx ' ' 



(35) 



— — = on X d^z , 6 — Q on dflx x il^ 
oz 



where g E L'^{fl) is a given function, with mean value along the z-coordinate equal to zero. 
g = 0. Moreover we shall explain why we have to introduce the Lagrange multiplier in order 
to solve numerically this equation. We remark that in contrast to section [3] we omitted 
for simplicity reasons the primes for 0, which indicated the fluctuation functions with zero 
mean value. 

The weak form of ([35]) is 

a((/.»=e(5,^), V^eV, (36) 

or equivalently 

m((/),V)=£(g,V), VV^eV, (37) 

where m{-,-) is the coercive bilinear form defined in (j25p . Let us now consider the corre- 
sponding discrete problem 

a[(j)h,iph) = eig.i^h) , VV'/iGV/i, (38) 

where the finite dimensional space Vh C V was introduced in section [XTl It can be seen that 
the property g = induces also in the discrete case that <f>h = Thus, following the same 
arguments as for the continuous case, we can show that equation psp is equivalent to 

m[cj)h,'iph) = £{g,i'h) , y-iph e Vh . (39) 

The Lax-Milgram theorem implies then the existence and uniqueness of a discrete solution 
(ph G V/i. The next theorem gives an estimate of the discretization error ||0 — (/);i||v. 
We shall suppose in the sequel, that the diffusion matrices A±, A^ and the function / are 
regular enough, to be able to use standard regularity/interpolation results. 
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Theorem 4.1 Let cj) eV be the unique solution of the continuous problem i36]) and (j)^ G Vh 
the unique solution of the discrete problem I138\) . Both solutions are elements of the normed 
space {U, II • \ \u), where 

U:^{ii-,-)eV with miu ■■= WdML-^m . 

Then we have the following discretization error estimate 

110 - 0/.|lv = \\dz<l> - dMh + ^\\d.<l> ~ d^Mlh < Ch^ , (40) 
with a constant C > independent of e > 0. Moreover, as (j), (ph G U, we have 

U-M\li<ch\ 

Proof: The fact that both solutions and (ph belong to the space U is an immediate 
consequence of the fact that the right-hand side of the equation (pH)) (rcsp. ([38])) satisfies 
g = 0. The discretization error estimate is rather standard. Denoting by 0/ the interpolant 
of (j) in the finite dimensional space Vh^ i-C. 

n=l k=l 

we have due to the coercivity of the bilinear form m{-, •) 

c||0 - 0/i|lv < '^('?^ - 0/1,0 - 0/i) = m{(l> - cj)h,(j>- pi) < c||(?!) - (^/illvll^ - 0/||v • 

Thus 

110 - 4>h\\v < c\\4> - PiWv ■ 

Standard Qi finite element interpolation results [26j yield for the interpolation error 

\\d.cj,~dMh + lia.0- 9.0/1 li. < c/i2(||a,,0||i. + p..0||i.) , 

and regularity results for the solution p of imply e^||9a;a;0||^2 + \ \dzz4>\\'^i^2 < ce^. This 
last estimate can be found by applying standard regularity results on the solution (p^ of 
the initial Singular Perturbation problem ([5|) (after a change of variable ^ := y/ex) and then 
exploiting the decomposition = p'^ + cpg. Thus, we have altogether with a constant c > 
independent of e > 

e||9x0-9.0/i||i2 + \\d,.p - d.^hWh < ch\ 



What is important to observe from the error estimate (|^n|) is that for e ^ the error 
110 ~ 0/i||//i in the standard £-independent _ff^-norm blows up. This is one argument why 
the Singular Perturbation model is inaccurate for £ <C 1. However, in the case where (p and 
(ph are elements of the space U, we have 1 10 — 0/i | |z^ < Ch^ independently of e, which means 
that we have convergence of the scheme in (U, \ \ ■ \ \u), uniformly in e > 0. The Poincare 
inequality implies then the uniform convergence in the || • ||l2 norm. The AP-scheme is thus 
equally accurate for every value of < £ < 1 . 

The discretization error <p ~ (ph is not the only error we are introducing when solving nu- 
merically ([55]) instead of ((551) . Indeed, ((55)) is nothing but a linear system 



Ala = V , 



(41) 
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to be solved to get the unknowns ank ■= 4>h{xn, Zk), where Vnk ■= £[g, Xni^k) and the discrete 
solution of (|38p is then reconstructed as 

(j)hix, -z) = X! X! ^rikXn{x)Kk{z) . 
n=l k=l 

Unfortunately the implementation of the system (|¥T|) introduces round-off as well as ap- 
proximation errors due for example to the numerical computation of a{xni^kT Xri^p 

). Thus 

the numerical resolution of (|4ip docs not yield the exact solution, but an approximation 
{ank)nk, solution of the slightly perturbed system 

Ma = w . (42) 

We are now interested in the error estimate — 0ft||v, as a function of the perturbation 
I — -5112, where || • ||2 denotes the Euclidean norm in ]R^=^^. 



Theorem 4.2 Let a be the exact solution of |^ j[ j and a the exact solution of the perturbed 
system Let 4>h S Vh and 4>h G Vh denote the corresponding functions 

<j>h{x,z) = ^^ankXn{x)K.k{z) , (ph{x,z) = '^'^ ankXn{x) Kk{z) . 
11=1 k=l n=l k=l 

Then we have 

e\\d,ci)h-dJh\\l2 + \\d,c^h~djH\\l2 < ^\\v ~ (43) 

with a constant c > independent of s > and h > 0. However, if both functions (j)h and 
belong to U, then we have the e-independent estimate 

Proof: Let us denote within this proof Enk ctnk — c^nk for n = 1, • • • ,-/Vx, k = 
1, • • • ,Nz and eh{x, z) := <ph{x, z) — (f>h{x, z), such that 



li^i ^) X] X] L^nkXn{x)Kk{z) ■ 



n=l k=l 

Moreover let := N^N^ and Y G be an arbitrary vector associated with the func- 
tion yh{x, z) = X)!^=i X^fcii ^nkXn{x)Kk{z). Thcn we have with (•, •)2 the cuclidcan scalar 
product in and M the discretization matrix of ((4T|) 

\\ME\\2 = sup — = sup — . 

Yes.'^ ,Y^o |K||2 yeR^.^T^o II^I|2 

Due to the fact that 

\\Y\\2<c\\yh\\L^<^\\yh\\v, 

we have 

I 7^1 1 m{yh,eh) /- m{yh,eh) _ 

\\ME\\2 = sup > cVe sup — — > cVe\\eh\\v ■ 

rem" ,Y^o iKlb ahev^.a^^o WVhWv 

Thus we get with a constant c > independent of e 

\\e,\\v<^\\ME\\2^^\\v-i}\\2. 



An asymptotic preserving scheme for strongly anisotropic elliptic problems 



20 



In the case the two functions (l)h and (j>h belong to U, i.e. Ch & U, we can exploit the fact 
that in Vl the Poincare inequality gives rise to \\Y\\2 < c||?/ft||L2 < c||y;i||iY. This yields, as 
m(-, •) is also coercive on U, that 

iisfrii m{yh,eh) , m{yh,eh) ,, ,, 
||M£'||2 = sup ' II > c sup — — >c\\eh\\u. 

yeB",Y5^0 Ik lb VhdU^yh^Q WUhWu 

and thus the e-independent estimate is proved. ■ 



Similarly as for the discretization error, we can deduce from the round-off error estimate 
([43|) that, for e ^ 0, the standard i?^-norm \ \4>h — </>/t||ffi explodes. However if we impose 
that both solutions (j)h and 4>h are elements of the space space of functions with mean 
value along the z-coordinate equal to zero, then we have the uniform estimate H^/t — |w ^ 
c\\v — v\\2 , and by the Poincare inequality \ \4)h — 0/i||l2 < c||w — v\\2. Unfortunately even if 
we know that (j)^ (zU, this is not necessarily true for 0^, if we discretize psp. But it can be 

achieved by forcing the numerical solution </ih, to satisfy (f>h = 0. Indeed, this can be done 
by introducing explicitly in the discrete problem (|38p the constraint iph = 0, such that it is 
much more ingenious to solve instead 

a{(t)h,'il'h) + b{Ph, rph) = sig, 4>h) , y-iph e Vh , 

(44) 

b{Qh,<l>h) = Q, yQheWh, 

where Wh C W was constructed in section 13.11 As mentioned in the continuous case this 
problem is equivalent for e > to the discrete problem ([55)1 . If (ph G Vh is the unique 
solution of dMl), then {(j)h,0) e x W;i is a solution of dH]). And if {(j)h,Ph) e V/i x W,, 
solves then = and (j)h € Vh is the unique solution of ([551) . This last statement 

is immediately proved by taking in the variational formulation (|44p only x-dependent test 
functions iphix) G Vh- By doing this, we can be sure that the numerical solution <j)h of ([44)) 

satisfies 4>h = 0, such that the error — (l)h\\u is uniformly bounded. This proves that 
the introduction of the constraint = in the AP-formulation is crucial and avoids the 
numerical difficulties associated with the original P-model. 



5 Conclusion 

In this paper we have introduced an Asymptotic Preserving formulation for the resolution of 
a highly anisotropic elliptic equation. We have shown the advantages of the AP-formulation 
as compared to the initial Singular Perturbation model and to its limit model, when the 
asymptotic parameter goes to zero. It came out that the AP-scheme is a powerful tool for 
the resolution of elliptic problems presenting huge anisotropics along one coordinate, and 
gives access to the simulation in a very easy and precise manner. The Asymptotic-Preserving 
method developed here relies on the decomposition of the solution in its mean part along the 
anisotropy direction, and a fiuctuation part. This integration along the anisotropy direction 
is easily performed in the context of Cartesian coordinate systems with one coordinate 
aligned with the direction of the anisotropy. In a forthcoming work [9] this procedure is 
extended to more general anisotropies. 
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